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1.0  SUMMARY 


The  presented  research  explores  a  new  approach  to  uncertainty  quantification  in  non-linear  systems 
with  application  to  orbit  trajectory  prediction  and  satellite  conjunction  analysis.  This  statistical 
approach  utilizes  novel  methods  to  build  better  uncertainty  state  characterization  in  the  context  of 
rare  event  prediction  while  keeping  computational  expense  low.  Current  models  are  not  well  suited 
to  predict  low  probability  regions  that  are  of  importance  to  accurately  calculating  conjunction  risk 
of  satellites.  Current  models  utilize  significant  but  sometimes  arbitrary  buffers  to  account  for  the 
unknown  true  statistical  distribution  of  satellite  position,  which  manifests  as  a  result  of  the  non¬ 
linear  properties  of  satellite  motion.  The  proposed  method  propagates  uncertainty  hypersurfaces, 
and  for  the  case  of  lower  dimensional  analysis  such  as  positional  uncertainty,  projects  these  hyper 
dimensional  surfaces  onto  lower  dimensions.  In  these  cases,  algorithms  are  utilized  to  find  the 
projected  uncertainty  contours  to  greatly  reduce  computation  costs.  These  algorithms  are  further 
optimized  by  exploiting  characteristics  of  the  non-linear  system  to  better  sample  the  uncertainty 
volume.  The  most  important  aspect  of  this  method  is  that  the  computation  cost  is  directly  related 
to  the  probability  and  desired  resolution,  making  it  computationally  efficient  for  characterization 
of  tail  probabilities.  This  method  can  be  used  to  assess  the  risk  of  rare  events  in  non-linear 
propagations,  and  is  here  applied  to  satellite  conjunction.  Furthermore,  this  method  will  be  applied 
to  the  inverse  problem  to  write  control  laws  for  tasking  ground  based  sensors  based  on  sensitivity 
analysis  of  uncertainty  quantification,  collective  multi  object  information,  and  collisional  risk 
assessments. 

2.0  INTRODUCTION 

In  recent  years,  the  number  of  tracked  objects  orbiting  earth  has  increased  significantly  due  to 
advancements  in  sensor  technology  and  processing  that  allow  for  smaller  and  smaller  objects  to 
be  monitored.  Currently,  over  40,000  objects  are  tracked  [45],  As  more  objects  are  detected  due  to 
increased  tracking  resolution,  the  risk  of  satellite  collision  is  better  predicted.  Active  satellites, 
decommissioned  satellites,  satellite  chunks,  and  other  space  junk  like  flakes  of  paint  and  lost  tools 
pose  a  threat  to  active  satellites,  as  seen  Figure  1 . 


Figure  1.  Satellite  Debris  Graphic 
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For  each  of  these  objects,  uncertainty  exists  in  location  and  velocity,  and  this  uncertainty  grows  as 
it  propagates  through  time  without  a  measurement  update.  When  two  propagated  uncertainty 
clouds  cross,  and  a  significant  conjunction  probability  is  associated  with  these  objects, 
surveillance  increases,  and  models  are  updated  to  estimate  collision  risk  [23-25],  This  risk 
estimation  is  not  very  precise,  and  the  error  associated  with  it  is  based  upon  linear  error  theory 
[1],  [47],  The  use  of  mean  position  or  maximum  likelihood  is  an  incomplete  picture  of  the 
uncertainty,  and  it  can  compromise  the  effectiveness  of  sensor  tasking  algorithms,  object 
identification,  and  calculation  of  collision  risk.  The  2009  collision  between  the  Iridium  33  and 
Cosmos  2251  satellites,  Figure  2,  shows  that  space  object  collision  is  a  real  threat  to  active 
satellites,  current  predictions  of  collision  risk  are  not  precise  enough,  and  sensor  tasking 
priorities  need  refining.  _ 


-1- 

T 

1 

t 

A 

1 

\  ■> 

w 

f-. 

Jl 

T 

-k 

Cu  i  mos  2251 

debris 

Figure  2.  Debris  Cloud  of  Colliding  Satellites 

Current  approaches  for  orbit  state  estimation  assume  the  probability  of  the  satellite  location  is 
Gaussian  in  3-dimensions,  and  propagates  it  as  growing  concentric  ellipsoids.  This  focus  on  esti¬ 
mating  the  maximum  likelihood  point  of  position  and  a  standard  deviation  measure  of  the  final 
distribution  is  important  for  identifying  objects  and  finding  them  after  several  orbits.  However, 
when  predicting  rare  events,  such  as  satellite  collisions,  rare  tail  probabilities  become  increasingly 
important.  Due  to  the  non-linearity  of  orbit  mechanics,  a  second  order  Gaussian  estimation 
encloses  a  lot  of  “dead  space”  where  the  satellite  has  a  zero  probability  of  being  located,  and  cuts 
out  regions  of  high  probability  density.  This  over-estimates  part  of  the  risk,  underestimates  another 
part,  and  does  not  keep  tabs  on  the  error  associated  with  the  risk  estimation.  By  improving  this 
model,  we  can  reduce  the  number  of  observations  required  to  tighten  the  uncertainty  volume,  can 
reduce  the  number  of  false-positive  conjunction  detections  which  often  require  costly  orbital 
maneuvers  that  limit  the  satellite’s  operational  lifetime,  and  can  reduce  false-negative  conjunction 
detections  which  result  in  a  collision,  mission  failure,  and  debris  clouds  that  are  difficult  to  track 
and  can  compromise  other  satellites. 
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This  work  tests  a  statistical  theory  that  may  more  accurately  determine  the  probability  density  of 
a  satellite’s  position  and  velocity  specifically  in  the  probability  tail  for  application  in  spacecraft 
conjunction,  debris  cloud  modeling,  and  other  non-linear  statistical  modeling  problems  [48].  A 
Monte  Carlo  scheme  is  used  as  a  “truth”  model  to  test  this  theory  when  analytical  methods  are  not 
available. 

Current  approaches  for  orbit  state  estimation  assume  the  probability  of  the  satellite  location  is 
Gaussian  in  3-dimensions,  and  propagates  it  linearly  [26-29].  Estimating  the  maximum  likelihood 
point  by  minimizing  variance  results  in  a  position  and  a  standard  deviation  measure  of  the  final 
distribution,  and  is  important  for  identifying  objects  and  finding  them  after  several  orbits. 
However,  when  predicting  events,  such  as  a  satellite  collision,  tail  probability  characterization 
becomes  increasingly  important.  Due  to  the  non-linearity  of  orbit  mechanics,  the  second  moment 
of  Gaussian  estimation  encloses  a  lot  of  “dead  space,”  where  the  satellite  has  a  zero  probability  of 
being  located,  and  cuts  out  regions  of  high  probability  density,  Figure  3.  The  figure  from  Reference 
[46],  demonstrates  that  initial  uncertainties  may  be  well  modeled  by  Gaussian  probability 
distributions,  but  quickly  diverge  from  these  elliptical  density  regions.  After  non-linear 
propagation,  second  order  Gaussian  estimation,  via  the  Extended  Kalman  Filter  (EKF)  and  the 
Unscented  Kalman  filter  (UKF),  will  enclose  a  lot  of  “dead  space”  where  the  satellite  has  a  zero 
probability  of  being  located,  and  cut  out  regions  of  high  probability  density.  This  overestimates 
part  of  the  risk,  underestimates  another  part,  and  does  not  keep  tabs  on  the  error  associated  with 
the  risk  estimation.  With  better  understanding  of  the  uncertainty  evolution  [13-15],  [49]  we  can 
construct  tracking  methodologies  that  reduce  the  number  of  observations  required  to  tighten  the 
uncertainty  volume,  can  reduce  the  number  of  false-positive  conjunction  detections  (which  often 
require  costly  orbital  maneuvers  that  limit  the  satellite’s  operational  lifetime),  and  can  reduce  false¬ 
negative  conjunction  detections  which  result  in  a  collision,  mission  failure,  and  debris  clouds  that 
are  difficult  to  track  and  that  can  compromise  other  satellites.  These  advanced  uncertainty 
quantification  techniques  can  then  be  used  to  efficiently  task  sensors  with  localization, 
disambiguation,  or  collisional  assessment  objectives  [3,4]. 
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Figure  3.  Using  Second  Order  Guassian  Models  on  Linear  and  Non-linear  Propagations  [5] 

The  process  of  collision  analysis  begins  after  observation  data  has  been  processed  to  form  orbit 
solutions  of  new  objects  and  maintain  solutions  of  preexisting  objects.  Typically,  satellites  are 
tracked  via  active  radar,  passive  radio  frequency  monitoring,  and  optical  tracking,  which  require 
image  processing  for  mean  and  covariance  measurements  as  well  as  satellite  identification.  [51], 
[52],  [53],  The  result  is  a  state  solution  (position  and  velocity)  that  typically  includes  an  associated 
covariance.  From  state  solutions,  mean  satellite  position  paths  are  propagated  to  identify  collision 
path  satellite  pairs  over  a  predefined  time  horizon.  At-risk  satellites’  state  pairs  are  reduced  to  a 
measure  of  mean  closest  approach  distance  at  a  particular  set  of  times.  While  this  determines  the 
most  likely  closest  approaches  it  may  ignore  greater  risks  that  involve  one  or  both  tail  probability 
overlap.  These  moments  of  closest  approach  generally  become  the  focus  of  collision  risk 
assessment.  Most  often,  collision  risk  is  assessed  by  assuming  the  velocity  is  negligible  and 
positions  are  Gaussian  post  propagation,  and  miss  distances  are  compared  to  the  added  diameters 
of  both  satellites  [28],  Instead  of  propagating  the  mean  points  and  adding  uncertainty  in  for  the 
moment  of  closest  approach,  this  research  is  interested  in  propagating  parts  or  all  of  the  uncertainty 
to  better  analyze  global  risk  of  collision  [31-37],  [39],  Current  methods  that  employ  this  technique 
exist,  such  as  Transformation  of  Variables  [2],  [7-12]  and  Monte  Carlo  methods  [2-5],  [28], 
[42-44],  but  are  computationally  expensive  [30],  Probability  Distribution  Function  (PDF)  evolution 
can  then  be  used  to  more  effectively  and  efficiently  task  ground  sensors  to  maintain  accurate 
location  data  for  satellites  approaching  collision  risk  parts  of  their  orbit,  and  to  best  close  large 
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uncertainly  estimates  as  well  as  increase  confidence  in  whether  or  not  satellite  operators  need 
to  maneuver.  Sensor  tasking  for  satellites  is  weighed  based  on  a  particular  satellite’s  required 
uncertainty  maximums,  the  expected  number  of  ground  stations  a  satellite  is  expected  to  pass 
during  its  orbit,  the  number  of  high  priority  satellites  also  crossing  a  particular  ground  station  in 
a  different  window  of  the  sky  during  the  same  time  interval,  and  the  expected  time  required  to 
find  the  satellite  [16  -22].  Energy  consumption  optimization  of  maneuvers  has  been  studied, 
however  it  does  not  account  for  future  collision  avoidance  [38].  Since  ground  stations  are  a 
limited  resource,  as  the  number  of  known  objects  in  orbit  increases  weekly;  there  is  a  need  for 
optimization  algorithms  of  sensor  measurements  and  tasking  priorities  [41]. 

3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

3.1  Analysis  of  Keplerian  Two-body  System  With  and  Without  Added 
Non-conservative  Forces 


Small  initial  uncertainty  distributions  in  velocity  and  position  of  satellites  quickly  grow  and  spread 
as  they  propagate  through  space.  The  characteristics  of  this  error  propagation  is  not  well 
understood.  A  breakdown  by  in-track,  cross-track,  and  radial  initial  error  uncertainties  may  help 
describe  the  satellite  position  and  velocity  probability  distribution  more  quickly  and  more 
accurately  than  current  models.  A  qualitative  analysis  was  done  on  these  propagations  and  are 
shown  and  described  in  this  section. 

To  analyze  the  characteristics  of  non-perturbed  satellites  with  position  and  velocity  uncertainty, 
we  used  a  Monte  Carlo  simulation  to  propagate  points  with  uncertainties  first  in  one  dimension, 
then  multiple  dimensions.  This  allowed  us  to  understand  the  fundamental  propagation 
characteristics  and  then  understand  how  they  mix.  Gaussian  noise  with  various  uncertainty 
characteristics  was  added  to  the  measured  nominal  orbit,  and  then  run  through  the  propagation 
program  to  calculate  the  position  and  velocity  vectors  of  each  Monte  Carlo  point  at  time,  t, 
relative  to  the  mean  orbit  period,  T  P. 

In  order  to  create  a  three  dimensional  probability  space  shape,  we  need  to  be  able  to  propagate 
Monte  Carlo  points  to  a  particular  time  in  their  non-perturbed,  Keplerian  orbit.  There  are  many 
methods  to  calculating  the  position  of  a  satellite  given  initial  conditions  including  Kepler’s 
problem,  Gooding’s  method,  and  Double-R.  Six  independent  pieces  of  information  are  required 
to  construct  an  orbit  [40,  50].  These  can  be  position  vectors,  velocity  vectors,  and  time  of 
observations.  For  simplicity,  we  programmed  a  simple  algorithm  that  inputs  a  single  radar 

observation  with  values  for  p,  p,  El,  E  l,  Az,  and  A  z,  (distance  from  radar  location  to  satellite, 
velocity  along  the  distance  vector,  elevation  angle,  change  in  elevation  angle,  azimuth  angle 
measured  clockwise  from  North,  and  change  in  azimuth  angle,  respectively),  and  combines  that 
with  the  radar  site  location  defined  with  6,  <p,  and  H,  (local  sidereal  time  or  angle  between  vernal 
equinox  and  radar  site  location,  latitude  angle,  and  altitude  of  radar  location  above  mean  sea  level, 
respectively),  calculates  the  position  and  velocity  vectors  in  IJK  coordinates,  and  then  uses 
Kepler’s  problem  to  solve  for  final  position  and  velocity  vectors  at  a  time,  t. 
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Noise  can  be  added  to  our  single  radar  observations  in  a  variety  of  different  ways,  i.e.  to  the 
South-East-Zenith  (SEZ)  position  and  velocity  vectors,  to  the  final  UK  position  and  velocity 
vectors,  etc.  In  order  to  separate  different  characteristics  of  the  final  distribution  shape,  we  added 
noise  to  the  initial  position  vector  with  respect  to  the  center  of  the  earth,  rO  and  the  initial 
velocity  vector  with  respect  to  the  center  of  the  earth,  vO,  in  the  following  axes:  in-track,  vO 
direction;  radial,  rO  direction;  and  cross-track,  vO  x  rO  direction.  The  velocity  standard  deviation 
was  added  as  a  portion  of  the  magnitude  of  the  initial  velocity,  vO,  and  the  position  standard 
deviation  was  added  as  a  portion  of  the  magnitude  of  the  position  vector,  rO.  These  uncertainty 
axes  are  not  necessarily  completely  independent,  since  the  velocity  and  position  vectors  are  not 
necessarily  perpendicular  to  each  other.  For  simplicity,  we  have  chosen  an  orbit  where  the 
position  and  velocity  are  near  perpendicular. 


3.1.1  Error  Propagation  in  the  In-Track  Direction 

For  velocity  uncertainty  in  the  in-track  direction  only:  points  begin  all  together  at  the  initial  position 
vector.  Points  with  higher  velocities  will  orbit  on  a  slower,  larger  ellipse,  while  points  with  lower 
velocities  will  orbit  on  a  quicker,  smaller  ellipse.  The  result  is  that  the  points  spread  out  along  the 
orbit  track,  the  slower  velocity  points  trailing  behind,  the  quicker  velocity  points  leading.  The 
slower  trailing  points  have  a  smaller  distance  to  travel  to  complete  one  revolution,  so  they  will  switch 
to  become  leading  points,  while  the  outer  points  which  have  a  larger  orbit  become  trailing  points. 
The  line  of  points  is  two-dimensional  in  the  plane  of  the  mean  orbit.  The  inner,  now  leading  points 
curl  in  to  the  initial  point,  while  the  now  trailing  outer  points  spread  out  along  their  larger  orbits. 
As  time  goes  on,  these  points  form  an  elliptic  swirl.  The  propagation  of  in-track  uncertainty  has 
the  greatest  impact  on  the  final  shape  of  the  total  probability  distribution.  Figure  4  and  Figure 
5  show  the  evolution  of  this  phenomena  with  and  without  position  noise.  Figure  4  demonstrates  a 
velocity  standard  deviation  at  5%  in  the  in-track  direction  and  zero  in  all  other  velocity  and  position 
directions.  The  blue  line  represents  the  initial  mean  position  vector,  and  the  red  line  represents  the 
final  mean  position  vector  with  propagation  in  the  clockwise  direction.  Figure  5  demonstrates  a 
velocity  standard  deviation  at  5%  in  the  in-track  direction  and  1%  in  all  position  directions.  The 
blue  line  represents  the  initial  mean  position  vector,  and  the  red  line  represents  the  final  mean 
position  vector,  propagation  is  in  the  clockwise  direction. 
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Figure  4.  In-Track  Velocity  Uncertainty 


Figure  5.  In-Track  Velocity  Uncertainty  with  Position  Uncertainty 
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3.1.2  Error  Propagation  in  the  Radial  Direction 


For  velocity  uncertainty  in  the  radial  direction  only:  points  begin  all  together  at  the  initial  position 
vector.  Points  with  positive  radial  velocity  uncertainty  will  complete  half  their  first  orbit  more 
slowly  than  points  with  negative  radial  velocity  uncertainty,  and  vice  versa  for  the  second  half  of 
the  first  orbit.  For  situations  where  the  initial  position  and  velocity  vectors  are  perpendicular,  or  if 
the  mean  orbit  is  circular,  the  mean  orbit  will  have  the  least  energy  and  complete  its  orbit  first. 
Points  with  positive  radial  velocity  uncertainty  will  trail  behind  the  mean  orbit  outside  of  it,  while 
points  with  negative  radial  velocity  uncertainty  will  lead  the  mean  orbit  within  it.  The  points  will 
then  switch  from  inside  to  outside  180  degrees  from  the  initial  point.  The  now  inner  points  that 
trailed  behind  will  start  to  catch  up  to  the  mean,  while  the  newly  outer  points  that  lead  the  way  will 
slow  down.  The  mean  point  will  always  complete  its  orbit  first  and  the  other  points  will  trail  more 
and  more  behind.  This  two-dimensional  line  within  the  plane  of  the  nominal  orbit  will  curl 
completely,  with  the  mean  point  leading  the  way  at  t  =  IT P,  the  orbital  period  for  the  nominal 
orbit.  Figure  6  and  Figure  7  show  the  evolution  of  this  phenomena  with  and  without  position  noise. 
Figure  6  demonstrates  a  velocity  standard  deviation  at  5%  in  the  radial  direction  and  1%  in  all 
position  directions.  The  blue  line  represents  the  initial  mean  position  vector,  and  the  red  line 
represents  the  final  mean  position  vector,  propagation  is  in  the  clockwise  direction.  Figure  7 
demonstrates  a  velocity  standard  deviation  at  5%  in  the  radial  direction  and  zero  in  all  other 
velocity  and  position  directions.  The  blue  line  represents  the  initial  mean  position  vector,  and  the 
red  line  represents  the  final  mean  position  vector,  propagation  is  in  the  clockwise  direction. 
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Figure  6.  Radial  Velocity  with  Position  Uncertainty 
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Figure  7.  Radial  Velocity  Uncertainty 
3.1.3  Error  Propagation  in  the  Cross-Track  Direction 

For  velocity  uncertainty  in  the  cross-track  direction  only:  points  begin  all  together  at  the  initial 
position  vector.  The  velocities  farther  away  from  the  mean  will  add  energy  to  the  orbit  and 
therefore,  will  complete  their  revolution  more  slowly.  The  mean  point  will  lead  the  way,  with  the 
other  points  trailing  behind  more  and  more.  Since  the  orbits  for  all  points  begin  at  the  same  point, 
the  orbital  plane  for  all  orbits  will  cross  at  the  initial  point  and  also  at  the  opposite  point.  However, 
these  crossings  will  occur  at  different  times  due  to  the  difference  in  orbital  period  for  points  with 
more  cross-track  uncertainty.  The  mean  point  will  cross  the  nodal  point  first,  followed  by  the  rest 
of  the  points  who  switch  from  one  side  of  the  cross-track  mean  orbit  to  the  other.  This  results  in  a 
little  loop  which  unfurls  to  a  curve,  then  comes  back  to  a  point  right  when  the  mean  orbit  reaches 
the  next  node,  and  all  the  points  switch  directions  in  succession.  This  three  dimensional  curve 
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roughly  follows  the  curvature  of  the  planet  for  circular  mean  orbits.  Figure  8  and  Figure  9  show 
the  evolution  of  this  phenomena  with  and  without  position  noise.  Figure  8  demonstrates  a  velocity 
standard  deviation  at  5%  in  the  cross-track  direction  and  1%  in  all  position  directions.  The  blue 
line  represents  the  initial  mean  position  vector,  and  the  red  line  represents  the  final  mean  position 
vector,  propagation  is  in  the  clockwise  direction.  Figure  9  demonstrates  a  velocity  standard 
deviation  at  5%  in  the  cross-track  direction  and  zero  in  all  other  velocity  and  position  directions. 
The  blue  line  represents  the  initial  mean  position  vector,  and  the  red  line  represents  the  final  mean 
position  vector,  propagation  is  in  the  clockwise  direction. 


Figure  8.  Cross-Track  Velocity  Uncertainty  with  Position  Uncertainty 


Approved  for  public  release;  distribution  is  unlimited. 
10 


Figure  9.  Cross-Track  Velocity  Uncertainty 
3.1.4  Error  Propagation  in  All  Directions 

For  velocity  uncertainty  in  all  directions:  Points  begin  all  together  at  the  initial  position  vector. 
The  points  then  spread  out  with  in-track,  cross-track,  and  radial  propagation  characteristics.  The 
points  begin  to  form  a  type  of  rope  wrapping  around  themselves  as  would  be  seen  with  the  in-track 
swirl.  They  form  several  nodal  points  whenever  they  cross  the  initial  position  vector,  or  its  opposite 
at  each  half  orbit  beginning  at  1.5  initial  orbits.  There  is  a  tilt  associated  with  each  anti-nodal  point 
depending  on  the  initial  cross-track  and  radial  standard  deviations,  giving  the  distribution  a  twisted 
ribbon-like  look.  Figure  10  and  Figure  11  show  the  evolution  of  this  phenomena  with  noise  and 
the  evolution  of  this  phenomena  without  noise  at  varying  angles.  Figure  10  demonstrates  a 
standard  deviation  at  5%  in  each  velocity  direction  and  1%  in  all  position  directions.  The  blue  line 
represents  propagation  is  in  the  clockwise  direction.  Figure  1 1  demonstrates  a  standard  deviation 
at  5%  in  each  velocity  direction  and  zero  in  all  position  directions.  These  graphs  are  taken  at  2.25. 
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TP  and  have  been  rotated  to  help  show  the  structure  of  the  distribution.  The  blue  line  represents 
the  initial  mean  position  vector,  and  the  red  line  represents  the  final  mean  position  vector.  The 
right  column  viewing  angle  is  looking  directly  down  the  end  of  the  initial  mean  position  vector, 
propagation  is  in  the  clockwise  direction. 


Figure  10.  Velocity  and  Position  Uncertainty 
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Figure  11.  Velocity  Uncertainty  at  Different  Angles 
3.1.5  Analysis  of  Keplerian  Motion 

The  objective  is  to  be  able  to  completely  describe  the  3-dimensional  probability  distribution  with 
the  least  amount  of  independent  information.  This  could  be  used  to  mathematically  determine  the 
probability  distribution  without  the  use  of  propagation,  as  well  as  to  connect  points  that  lie  on 
particular  cumulative  probability  distribution  surfaces.  Even  without  completely  determining  the 
probability  shape  and  structure,  this  uncertainty  knowledge  could  be  used  to  better  estimate  the 
probability  clouds  of  satellites  with  uncertainties.  This  could  result  in  a  more  accurate 
approximation  to  satellite  position  probabilities  which  be  valuable  for  cursory  and  initial 
conjunction  analysis. 

Nine  properties  of  the  probability  structure  were  determined  from  the  above  Monte  Carlo 
simulations.  These  properties  may  or  may  not  be  independent  and  may  or  may  not  be  able  to 
describe  the  full  3-dimensional  probability  distribution. 
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In-Track  velocity  uncertainty  will  affect: 

•  In-plane  orientation 

•  Curvature 

•  Length 

Radial  velocity  uncertainty  will  affect: 

•  Radial  width 

•  Radial  substructure 

•  Radial  orientation 

Cross-Track  velocity  uncertainty  will  affect: 

•  Cross-Track  width 

•  Cross-track  tilt 

•  Cross-track  substructure,  substructure  of  tail 

This  analysis  gives  insight  to  how  uncertainty  propagates  through  time  for  satellites.  The  structure 
is  very  non-linear  which  requires  non-Gaussian  assessments  to  accurately  determine  the 
probability  of  satellite  location.  The  results  shown  here  describe  various  characteristics  of  the 
satellite  probability  structure  including  the  propagation  of  single  axis  uncertainties  as  well  as  the 
propagation  of  mixed  axes  uncertainties.  The  nine  characteristics  of  the  propagated  uncertainty 
described  here  could  fully  define  the  probability  structure  of  the  propagated  uncertainty  cloud. 
Since  the  initial  in-track  velocity  uncertainty  has  the  greatest  impact  on  the  final  shape  of  the 
probability  cloud,  a  preliminary  equation  has  been  written  to  describe  the  general  structure  of  in¬ 
track  uncertainty  propagation.  The  characteristics  along  with  equations  describing  this  motion 
could  be  used  to  quickly  and  accurately  assess  the  preliminary  probability  of  a  satellite’s  position 
at  any  point  in  time,  and  be  used  to  quickly  determine  potential  conjunction  situations.  Future  work 
includes  analyzing  these  characteristics  with  non-conservative  forces  added.  This  could  lead  to 
accurately  determining  the  satellite’s  position  probability  and  quickly  calculating  precise 
conjunction  probabilities. 

3.2  Studies  Conducted  for  Qualitative  Sensitivity  Analysis 


Orbital  uncertainty  and  debris  clouds  often  evolve  in  unexpected  ways  depending  on  the  initial 
conditions  and  perturbative  forces  acting  on  the  points.  An  understanding  of  this  is  critical  for 
optimizing  point  selection  and  accurately  applying  our  proposed  method  to  this  non-linear 
scenario.  For  instance,  points  with  initial  velocity  uncertainty  in  the  radial  direction  end  up  wholly 
contributing  to  uncertainty  in  the  in-track  direction  within  an  orbit.  Initial  position  uncertainty  in 
radial  and  in-track  directions  form  an  anti-node  that  flips  the  position  of  all  points  as  they  funnel 
through  it.  Without  a  careful  study  of  these  phenomena,  point  by  point  calculations  would  not  be 
possible. 
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When  the  only  uncertainty  is  position  or  velocity,  as  is  the  case  in  debris  fields  when  secondary 
collisions  (collisions  between  pieces  of  debris  after  the  initial  breakup)  are  assumed  negligible  and 
drag  forces  are  also  assumed  negligible,  this  method  would  seem  to  work  perfectly,  since  the  initial 
position  uncertainty  is  zero,  no  mixing  is  required  and  points  that  delimit  the  initial  velocity 
uncertainty  threshold  of  interest  remain  on  that  surface  throughout  propagation,  seen  in  Figure  12. 
Figure  12  plots  colors  points  that  are  within  1,  2,  3,  and  greater  than  3  standard  deviations  in 
position.  There  is  no  velocity  uncertainty  or  drag.  The  delimitations  between  the  colors  outlines 
standard  deviation  potentials.  There  is  no  color  mixing  when  there  is  only  error  in  one  type  of 
dimension,  in  this  case  position  uncertainty. 
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Figure  12.  Shape  Evolution  with  Colored  Standard  Deviation  in  Position 
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Understanding  motion  under  two-body  motion  is  helpful  for  optimizing  initial  point  selection 
for  use  with  this  method  and  for  understanding  how  initial  uncertainties  affect  the  final  shape  of 
the  uncertainty  distribution.  Figure  13  and  Figure  14  show  how  different  components  of 
uncertainty  affect  the  propagation  shape  of  a  satellite  for  a  mean  orbit  that  is  slightly  elliptical. 
Figure  13  shows  the  shape  evolution  from  t  =  0  to  t  =  ,5T P  in  AT P  increments,  where  T P  is  the 
nominal  period.  5000  Monte  Carlo  points  following  2-dimensional  2-body  motion  with  no 
perturbations  were  assigned  initial  uncertainties  as  labeled  where  5%  of  the  initial  position 
vector  magnitude  was  added  for  position  uncertainty  and  5%  of  the  initial  velocity  vector 
magnitude  was  added  for  velocity  uncertainty.  Interesting  phenomena  can  be  seen  in  unexpected 
curls,  rotations  and  pinch  points.  Figure  14  shows  the  shape  evolution  from  t  =  0  to  t  =  ,5T  P 
in  AT  P  increments,  where  T  P  is  the  nominal  period.  5000  Monte  Carlo  points  following  2- 
dimensional  Two-body  motion  with  atmospheric  drag  perturbations  were  assigned  initial 
uncertainties  as  labeled  where  5%  of  the  initial  position  vector  magnitude  was  added  for 
position  uncertainty  and  5%  of  the  initial  velocity  vector  magnitude  was  added  for  velocity 
uncertainty.  Atmospheric  drag  was  highly  exaggerated;  atmospheric  density  was  assigned  as  8 
orders  of  magnitude  greater  than  actual  values.  While  the  central  position  of  the  distribution  has 
been  skewed  towards  Earth,  the  general  shapes  remain  intact. 
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Figure  13.  Shape  Evolution,  No  Perturbations 
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Figure  14.  Shape  Evolution,  with  Drag  Perturbations 

Uncertainty  in  radial  velocity  can  be  seen  to  contribute  primarily  to  the  in-track  uncertainty  after 
about  a  quarter  of  an  orbit.  This  is  important  in  picking  final  points  because  there  is  a  directionality 
that  needs  to  be  addressed  under  Keplerian  motion.  Since  we  are  combining  two  distributions,  a 
position  only  final  distribution  and  a  velocity  only  distribution,  with  the  option  to  do  this  with 
respect  to  a  mixed  distribution,  we  need  to  combine  points  that  are  in  the  same  direction  post 
propagation.  For  linear  propagations,  initial  and  final  points  in  the  same  direction  end  up 
contributing  in  the  same  direction  post  propagation.  This  is  not  true  of  Keplerian  motion.  Initial 
velocity  points  need  to  be  chosen  90°  offset  from  position  uncertainty  in  order  to  properly  combine 
with  position  uncertainties  for  propagations  longer  than  a  quarter  of  an  orbit.  This  can  be  seen 
when  looking  at  initial  radial  velocity  uncertainty  and  in-track  position  uncertainty.  All  the  points 
line  up  best  in  this  combination  meaning  that  the  uncertainties  are  contributing  to  the  same  final 
phenomenon. 
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Other  interesting  features  that  can  be  seen  are  a  flip  in  position  uncertainties  at  about  .3T  P  .  The 
resulting  node  is  best  seen  when  velocity  uncertainty  is  not  included.  This  flip  corresponds  to  the 
outside  points  with  greater  initial  potential  energy  converting  that  potential  energy  to  kinetic 
energy  and  catching  up  to  the  inside  points  with  a  faster  velocity.  This  flip  will  undo  itself  on  the 
other  side  of  the  orbit  and  will  continue  this  oscillation  if  no  other  forces  are  present.  When  only 
velocity  uncertainty  is  included,  a  single  node  appears  at  the  origin  of  the  points.  All  points  will 
filter  through  this  node  as  they  complete  their  orbits. 

3.3  Study  on  Post-propagation  Covariance  Calculations 


Moment  estimation,  while  mathematically  easy  to  calculate,  does  not  completely  describe  non- 
Gaussian  distributions;  two  distributions  might  have  identical  moments,  but  have  very  different 
distributions.  For  very  non-linear  distributions,  moments  estimations  become  less  indicative  of  the 
overall  distribution  and  can  be  an  extremely  poor  representation  of  the  distribution,  see  Figure  15. 
Figure  15  shows  the  statistical  mean  and  standard  deviation  were  calculated,  along  with  UKF 
approximation,  and  the  statistical  mean  and  standard  deviation  of  points  on  the  basis  loops 
projected  via  the  proposed  method.  None  of  these  moment  estimations  do  a  very  good  job  of 
describing  the  distribution.  The  mean  points  do  not  match  up  nor  do  they  indicate  a  point  of 
significant  probability.  Moments,  while  useful  in  some  applications,  are  not  a  helpful  way  of 
describing  non-linear  distributions. 


8000  r 


6000  - 


Keplerian  Motion  UKF  Calculation 


£ 

u 
c 

EU 

-2000 

~o 

-4000  k 


-10000 


Monte  Carlo  points 
X  Monte  Carlo  mean 
—  Monte  Carlo  Covariance 


Basis  Loops 
Basis  Loop  mean 
- Basis  Loop  Covariance 

2  UKF  sigma  points 
UKF  mean 
UKF  Covariance 


-0,5  0 

x  distance  (km) 


xio* 


Figure  15.  Moment  Estimation  Comparison 
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For  nonlinear  distributions,  we  are  interested  in  a  description  of  the  distribution  that  is  useful  for 
calculating  our  satellite  tracking  objectives  (collision  risk  assessment,  localization, 
disambiguation,  etc.).  We  want  a  quantitative  description  of  the  uncertainty,  in  our  case  a  set  of 
final  uncertainty  loops,  calculated  in  such  a  way  that  each  uncertainty  loop  encompasses  a 
particular  proportion  of  the  distribution.  Furthermore,  any  uncertainty  loop  calculated  in  this  way 
must  be  completely  contained  by  all  uncertainty  loops  (also  calculated  by  this  method)  that 
encompass  a  greater  proportion  of  the  distribution  (i.e.  the  loops  are  concentric,  are  larger  if  the 
uncertainty  enclosed  is  larger,  and  do  not  cross  each  other).  These  restrictions  provide  a 
mathematical  fidelity  for  future  calculations. 

For  propagations  of  distributions  with  initial  Gaussian  uncertainties,  we  have  found  that  the 
propagation  of  the  full  dimensional  uncertainty  hyper  ellipse  through  non-linear  functions  meets 
these  mathematical  objectives  post  propagation  [6].  The  projection  of  the 
multidimensional  uncertainty  loop  onto  lower  dimensional  planes  meets  these  mathematical 
conditions  for  quasi-linear  propagations,  and  approximates  them  for  highly  non-linear 
propagations.  There  are  at  least  four  ways  to  estimate  these  projections: 

1.  Full  projection 

2.  Basis  loop  projection 

3.  Dimension  reduced  calculations 

4.  Edge  loop  calculation 

3.3.1  Full  Projection 

Beginning  with  n-dimensional  Gaussian  uncertainty,  we  sample  points  as  densely  as  required  (or 
limited  by  computation  resources)  on  the  surface  of  the  n-dimensional  ellipsoid  that  denotes  a 
particular  uncertainty  threshold.  These  points  are  sent  through  the  non-linear  propagator  and  then 
projected  onto  an  m-dimensional  surface  (where  m<=n ).  In  the  case  of  2-D  Keplerian  motion, 
we  have  four  initial  uncertainty  dimensions  (x,  y,  x ,  and  y-),  and  we  are  projecting  onto  two 
dimensional  surface  (for  instance  x  and  y  position,  or  x  and  y  velocity,  etc.)  Unfortunately  this 

requires  sn  number  of  points,  where  s  is  the  number  of  sample  points  per  dimension.  This  number, 
s,  can  be  as  low  as  2,  however  this  quickly  becomes  unfeasible  for  high  dimensional  systems  and 
can  be  comparable  to  the  number  of  samples  required  or  Monte  Carlo  simulations  (The  general 

rule  of  thumb  for  Monte  Carlo  sample  size  is  1 0n).  If  the  computational  resources  are  available, 
this  method  gives  high  fidelity  estimates  of  the  uncertainty  loop  even  at  extremely  low  thresholds. 
While  a  Monte  Carlo  simulation  will  have  lower  accuracy  and  require  many  more  points  to  reliably 
calculate  tail  probabilities,  this  method  will  give  accurate  calculations  even  at  tail  probabilities.  The 
benefits  of  this  method  also  allow  the  sampling  of  additional  independent  points,  r,  in  regions  that 

require  additional  resolution  (beyond  the  resolution  given  by  sn  sample  points)  without  the  need  to 
resample  the  entire  distribution.  In  addition,  this  method  allows  accurate  calculations  of  the  error, 
related  to  the  local  resolution,  and  is  directly  proportional  to  n. 
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For  Monte  Carlo  tail  probability  calculations: 


Number  of  samples  ~  (10  +  p)  v  J 

Where,  p  is  a  number  which  is  sufficiently  large  to  allow  a  high  probability  of  accuracy  in  the  tail 
region. 

Full  projection  tail  probability  estimation: 

(2) 

Number  of  samples  ~  s  +  £ 

Where  5  is  a  number  greater  than  1  (recommended  2-5  to  start),  and  £  is  the  length  of  the  perimeter 
of  the  projected  hypersurface  where  each  dimension  is  scaled  by  its  resolution  requirement. 

This  method  is  good  when  the  system  is  completely  unknown,  however  it  requires  a  large  number 
of  extra  points  to  be  projected  before  resolutional  points  can  be  applied.  It  is  a  significant 
improvement  to  Monte  Carlo  computational  requirements,  however  can  get  out  of  hand  for  high 
dimensions,  or  for  long  highly  non-linear  propagations  that  go  through  many  stretching/folding 
episodes. 


3.3.2  Basis  Loop  Projections 

Another  way  to  sample  is  by  the  use  of  Basis  loops.  These  loops  are  the  intersection  of  the  full 
dimensional  hyper  ellipse  to  sets  of  two-dimensional  orthogonal  planes  each  of  which  includes  the 
mean  point.  For  our  purposes,  we  simply  used  planes  formed  by  each  uncertainty  dimension 
combination.  Figure  16  shows  a  4-dimensional  hyper  ellipse,  and  its  six  2-D  Basis  loops  projected 
into  three  dimensions.  Figure  16  shows  sample  uncertainty  for  an  arbitrary  Keplerian  2-D  motion 
problem.  The  basis  loops  are  as  follows:  px  vs  py  in  red,  vx  vs  vy  in  green,  px  vs  vx  in  marigold, 
py  vs  vy  in  blue,  px  vs  vy  in  orange,  and  vx  vs  py  in  turquoise,  where  p  denotes  position  and  v 
denotes  velocity.  Some  of  the  basis  loops  were  projected  down  to  lines  in  the  px-py  -vx  dimensional 
space  because  as  ellipses,  their  other  uncertainty  dimension  was  in  vy. 

For  this  4-D  example,  we  could  also  construct  four  1-D  basis  loops  (lines),  six  2-D  loops  (ellipses), 
or  four  3-D  basis  loops  (spheres).  The  six  2-D  basis  loops  give  the  most  information  per  number 
of  propagated  points,  i.e.  less  points  would  be  propagated  for  the  six  1-D  lines,  but  the  projection 
outline  would  be  poorly  defined,  and  while  the  four  3-D  spheres  would  give  better  outline 
characteristics,  the  number  of  points  propagated  significantly  increases. 
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Initial  uncertainty  for  x  and  y  position  and  x  velocity 


Figure  16.  Initial  4-D  Hyper  Ellipse  with  Six  2-D  Basis  Loops  Projected  onto  3  Dimensions 

To  illustrate  this,  we  did  a  simple  2-D  polar  to  Cartesian  transformation.  Uncertainty  was  assigned 
in  all  four  dimensions  of  polar  uncertainty,  r,  9,  r',  and  O' ,  basis  loops  were  constructed,  and  then 
transformed  to  Cartesian  coordinates  via  the  following  transformation  function: 


X 

r  cos  6 

y 

r  sini9 

X 

r  cos  $  —  r$  sin  9 

.  y  \ 

r  sin  $  +  rO  cos  $ 

The  resulting  basis  loops  were  then  plotted  in  Cartesian  coordinates,  as  seen  in  Figure  17. 
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(a)  2-D  Basis  Loop  Projection 
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(b)  3-D  Basis  Loop  Projection 


Figure  17. 2-D  Polar  to  Cartesian  Transformation  Projections  onto  2-D  Orthogonal  Planes 

Figure  17  shows  Uncertainty  in  r,  9,  f ,  and  O'  were  transformed  to  Cartesian  coordinates,  then 
projected  onto  two  dimensional  orthogonal  planes:  (a)  shows  all  six  2-D  basis  loops  projected  onto 
each  axes  pair,  and  (b)  shows  these  same  basis  pairs  along  with  a  the  four  combinations  of  3- 
dimensional  spheres,  propagated  and  projected  onto  each  two  dimensional  axes  pair.  The  3-D  basis 
loops  fill  out  closer  to  the  edges  of  the  actual  probability  outline  than  the  2-D  basis  loops,  but 
require  many  more  points  to  do  so.  Ideally,  propagation  would  be  limited  to  only  points  that  end 
up  lying  on  this  outline,  without  propagating  all  the  superfluous  inside  points. 
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3.3.3  Edge  Loop  Propagation 


Only  points  that  lie  on  the  boundary  of  the  propagated  and  projected  (to  m-dimensions)  hyper 
ellipse  are  strictly  necessary.  With  a  knowledge  of  the  propagation  function,  only  the  points  that 
lie  on  an  m-dimensional  loop  that  share  the  edge  with  the  n-dimensional  full  projection  would  be 
propagated.  This  would  limit  the  sample  size  to  sm.  This  edge  loop  may  not  be  continuous  in  m- 
dimensional  space,  as  folds  in  the  distribution  can  produce  non-convex  shapes. 

The  sensitivity  analysis  that  is  ongoing  will  allow  dimension  reduction  in  calculations,  and 
possibly  allow  for  the  smart  sampling  of  this  piecewise  edge  loop  with-  out  projecting  the  entire 
n-  dimensional  surface.  This  would  need  to  be  analyzed  for  each  non-linear  system  and  will  be 
analyzed  for  Keplerian  motion  in  Section  3.4.  Furthermore,  with  knowledge  of  the  final  function 
sets  and  the  edge  loop  points,  projection  of  the  minimum  number  of  points  required  to  fully  de¬ 
scribe  a  function  in  that  set  could  greatly  reduce  the  number  of  sample  points  required.  For 
instance,  three  points  are  required  to  fully  describe  a  circle  or  a  parabola,  given  any  three 
independent  points  on  these,  we  can  fully  describe  the  final  function. 

3.3.4  Edge  Loop  Solver 

While  the  edge  loop  propagation  requires  high  knowledge  of  the  propagation,  an  edge  loop  solver 
would  need  little  to  no  prior  knowledge  of  the  system  to  run.  An  edge  loop  solver  could  start 
with  low  dimensional  basis  loops,  and  then  solve  for  extrema  points  on  the  sections  of  the  basis 
loops  that  form  the  outer  edge,  i.e.  those  sections  that  are  not  contained  within  any  other  loop.  The 
initial  positions  of  these  points  would  then  be  connected  along  great  circles  of  the  initial  hyper 
ellipse,  and  these  arcs  would  be  sampled,  propagated,  and  projected.  Extrema  from  these  arcs  could 
also  be  connected  by  initial  great  circles  and  propagated  until  the  change  in  border  points  is  below 
the  resolution  requirements.  Alternately,  points  near  these  great  circle  arcs,  with  some  cross¬ 
dimensional  term,  could  also  be  propagated  to  see  if  they  expand  the  edge  boundary  or  not. 
Algorithms  built  in  these  or  similar  ways  could  solve  for  the  boundary  of  the  propagated  and 
projected  full-dimensional  hyper  ellipse.  While  this  may  require  several  iterations,  the  number  of 
sample  points  would  still  remain  low,  and  various  checks  could  be  implemented  to  determine 
accuracy.  Edge  loop  solvers  could  also  be  seeded  with  initial  loops  that  are  known  to  the  system 
for  producing  high  uncertainties.  Some  uncertainties  or  combinations  of  uncertainties  could  lead 
to  minimal  impact  on  the  final  distribution,  and  if  these  are  found,  they  can  help  sample 
intelligently  for  decreased  computational  burden. 

3.4  Basis  Loop  Examples  for  Keplerian  Motion  in  2-D 

All  basis  loops  were  given  initial  uncertainty  for  in-track  and  radial  position  and  velocity.  Figure 
18  shows  how  basis  loops  form  sets  that  span  the  uncertainty  space.  These  pairs  can  give  insight 
into  how  the  uncertainty  behaves  as  well  as  where  to  look  for  the  actual  full  uncertainty  border.  In 
Figure  18,  basis  loops  form  sets  that  span  the  full  dimensional  uncertainty.  In  this  case  of  a  simple 
spring,  px  vs  py  in  red  and  vx  vs  vy  in  green  form  a  set,  px  vs  vx  in  marigold  and  py  vs  vy  in 
blue  form  a  set,  and  finally,  px  vs  vy  in  orange  and  vx  vs  py  in  turquoise  form  a  set,  where  p  denotes 
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position  and  v  denotes  velocity.  The  magenta  loop  shows  the  actual  uncertainty  loop.  This  loop 
can  be  found  by  any  set  of  base  pairs  by  summing  the  variances  as  is  expected  for  independent, 
identically  distributed  random  variables.  Loops  that  are  flattened  (and  therefore  bound  zero  area) 
trace  a  transformed  coordinate  axis.  The  direction  of  the  collapse  indicates  information  that  does 
not  influence  the  final  distribution.  The  corresponding  loop  pair  will  then  extend  to  the  actual 
uncertainty  in  that  transformed  direction.  The  smaller  an  uncertainty  loop  by  length,  the  less 
influence  it  has  on  the  final  uncertainty,  and  the  closer  its  pair  is  to  the  actual  uncertainty. 

Figure  19  shows  an  example  of  how  basis  loops  propagate  in  Keplerian  two-body  motion,  and  how 
we  can  use  knowledge  of  these  to  understand  the  underlying  dynamics.  In  Figure  19,  the  basis 
loops  are  as  follows:  px  vs  py  in  red,  vx  vs  vy  in  green,  px  vs  vx  in  marigold,  py  vs  vy  in  blue,  px 
vs  vy  in  orange,  and  vx  vs  py  in  turquoise,  where  p  denotes  position  and  v  denotes  velocity.  The 
red  loop,  (paired  with  the  green  loop)  is  almost  completely  collapsed,  indicating  that  it  does  not 
influence  the  final  uncertainty  along  the  non-linear  axis  of  collapse.  Along  this  axis,  the  green  loop 
stretches  out  to  meet  the  actual  uncertainty.  This  pair  indicates  an  uncertainty  combinations  that 
add  together  for  great  influence  one  way,  and  that  work  against  each  other  to  reduce  influence  in 
another  way.  This  might  be  helpful  to  collapse  the  problem  to  a  lower  dimensional  space  without 
compromising  the  estimation  integrity. 

Figure  20  shows  how  an  edge  loop  solver  might  sample  points.  In  Figure  20,  the  basis  loops  are 
as  follows:  px  vs  py  in  red,  vx  vs  vy  in  green,  px  vs  vx  in  marigold,  py  vs  vy  in  blue,  px  vs  vy  in 
orange,  and  vx  vs  py  in  turquoise,  where  p  denotes  position  and  v  denotes  velocity.  We  can  see  the 
loops  that  “stick  out”  of  the  rest  of  the  uncertainty  mass.  The  initial  positions  of  the  extrema  of 
these  protrusions  can  be  connected  by  great  circles  on  the  original  hyper  ellipse,  and  these  arcs 
sampled,  propagated,  and  added  to  the  uncertainty  estimate.  Extrema  of  these  arcs  can 
subsequently  be  reseeded  into  the  algorithm  to  push  the  boundaries  of  the  uncertainty  closer  to 
their  actual  values.  Tracing  the  calculated  boundaries  on  the  initial  hyper  ellipse  can  help  future 
sampling  by  indicating  which  uncertainty  combinations  are  most  important  to  propagate. 
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Figure  18.  Basis  Pairs 
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Figure  19.  Keplerian  Motion.  2-D  Basis  Loop  Projection  Early  in  Orbit 
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Figure  20.  Keplerian  Motion.  2-D  Basis  Loop  Projection  Late  in  Orbit 


3.5  Analytical  Equation  Propagation  for  Non-perturbed  Velocity 
Uncertainty 


The  shapes  shown  in  Figure  13  and  Figure  14  are  a  product  of  differential  equations  for  which  we 
do  not  have  analytical  solutions.  The  greatest  final  range  of  uncertainty  can  be  seen  as  originating 
from  in-track  velocity  uncertainty.  Describing  the  shape  of  this  curve  could  be  helpful  in 
circumventing  filtering  techniques  by  simply  jumping  to  the  final  solution.  However,  with  a 
numerical  chart  relating  eccentricity  to  true  anomaly,  the  error  can  be  minimized.  Figure  21  shows 
how  the  equation  derived  below  compares  to  the  actual  uncertainty  shape  using  a  second  order 
approximation  for  the  true  anomaly.  In  Figure  21  the  black  line  shows  the  actual  uncertainty  cloud 
after  10  orbits  for  initial  in-track  velocity  uncertainty  that  spans  from  95%  to  110%  of  the  mean 
velocity  magnitude  that  has  a  circular  orbit.  The  pink  line  shows  the  equation  result  using  a  second 
order  approximation  for  the  true  anomaly  in  time.  The  blue  line  represents  the  outline  of  the  Earth. 
Given  enough  velocity  uncertainty,  some  points  will  crash  into  the  earth. 
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Figure  21.  Comparison  of  Approximation  and  Exact  Solution 

Given  the  initial  nominal  measurement  vectors  for  position  and  velocity,  rO  and  vO,  the  positive 
and  negative  velocity  uncertainty,  vd+  and  v<5-,  and  some  way  to  identify  the  final  time  of  interest, 
tf.  Our  equation  explicitly  makes  use  of  tf,  but  any  of  the  following  can  be  converted  to  tf. 

1.  The  number  of  nominal  orbits  to  the  final  time  of  interest,  ntf 


tf 

TP0 


Tltf  -TPo 


27T  3 

:G0 


a0  = 


AH’ o 


vgr0 


2^t 


(4) 


Where  T  PO  is  the  nominal  orbital  period,  /u  is  the  gravitational  parameter  and  aO  is  the 
nominal  semi-major  axis. 
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2.  The  final  nominal  position  vector,  r/ 


COS  l/f 


eo  » 

e0  r/ 


(5) 


Where  eO  is  the  eccentricity  vector  of  the  nominal  point,  and  vf  is  the  angle  between  the 
final  nominal  position  and  nominal  periapsis.  This  vf  value  can  be  used  in  the  equations 
below  to  solve  for  tf . 

3.  The  angle  between  the  final  nominal  position  and  perigee,  vf 


tf 


[2&7T  +  (E(i/y)  -  e0sinE (vf))  -  (E(v0)  -  e0  sinE(i/0))] 


cosE(i') 


eQ  4-  cos  v 
1  +  e0  cos  v 


eo  • 

cost'o  = - 

eoro 


r0  -  (r0  •  v0)  v0 


(6) 


Where  k  is  the  number  of  times  the  nominal  position  passes  through  perigee  from  the 
initial  time  to  the  final  time  of  interest,  E(v)  is  the  eccentric  anomaly,  and  vO  is  the 
nominal  true  anomaly  at  initial  time. 

We  will  write  our  equation  as  a  function  of  the  magnitude  of  the  velocity,  v.  This  velocity  will 
range  from  vO  -  vd~  to  vO  +  vd+  since  the  uncertainty  is  all  in  the  velocity  direction,  the  uncertainty 
vector  is  in  the  same  direction  as  the  initial  velocity  vector.  We  assume  two-body  motion  with  no 
perturbations,  and  that  the  velocity  range  is  limited  to  motion  that  results  in  elliptical  orbits.  Each 
velocity  in  the  range  will  follow  its  own  elliptical  orbit  which  can  be  written  as: 


rM 


a  • (1  —  e2) 
1-1-6’  cos  {v) 


(7) 


Approved  for  public  release;  distribution  is  unlimited. 
29 


Where  r  is  the  radial  component  of  the  position  vector,  v  is  the  angular  component  of  the  position 
vector,  a  is  the  semi-major  axis,  and  e  is  the  eccentricity.  When  looking  at  the  collection  of  points 
with  different  initial  velocities,  all  points  will  follow  their  own  orbital  path  and  at  a  particular  point 
in  time,  tf  we  can  write  their  positions  as: 

Fbr  v0  -  v*_  <  v  <  v0  +  vs+ 

a(v)  -  (1  -  e(v)2)  (8) 

”V>  1  +e(v)  •  cos(t'(v,t/)) 


Where  r(v,  tf)  is  the  function  of  the  radial  di  stance  from  the  center  of  the  earth  to  points  with 
initial  velocity,  v  at  time  tf ,  a(v)  is  the  semi-major  axis  as  function  of  v,  e(y  )  is  the  eccentricity  as 
a  function  of  v,  and  v(v,  tf  )  is  the  angle  of  the  v  point  at  the  final  time,  tf .  We  now  have  to  find 
an  equation  for  a(y),  eiy)  and  v(v,  tf  ).  We  can  write  the  semi-major  axis,  a  in  terms  of  the  specific 
mechanical  energy,  E : 


2S 


Since  the  function  for  E  is: 


(9) 

(10) 


We  can  combine  Equations  9  and  10  in  terms  of  the  velocity,  v  and  the  initial  radial  magnitude, 
r0: 


a(v)  = 


(11) 


Which  can  be  rewritten  as: 


v2r0  —  2  fi 


(12) 


We  can  find  the  equation  for  e(v)  by  applying  the  following  initial  condition  to  Equation  4: 

When  z/(v,  tf)  —  r(v,  tf)=i0  (i3) 
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We  can  say  this  because  we  are  assuming  non-dispersional  forces,  so  all  orbits  will  cycle  back 
through  the  initial  radial  position  after  every  orbit.  We  can  apply  this  condition  to  Equation  8: 

a(v)  (l  -  e(v)2) 

rn  —  v  -  - \J—L  (14) 

1  +  e(v)  cos  {v  o) 

Can  be  rearranged  to: 

a(v)e(v)2  4-  r0cos(i/0)e(v)  +  r0  -  a(v)  =  0 


We  can  then  solve  for  e(v)  by  using  the  quadratic  formula: 


€ 


-r0cos(i>0)  ±  \/tqCos2{ 

>0 

)  4a(v)  (r0  -  a(v)) 

2a(v) 

(15) 


Since  the  eccentricity  is  less  than  1  for  an  ellipse,  we  can  make  the  following  constraint: 

When  cos  (i/0)  >  0, 

-r0cos(i>o)  +  \/rlcos2(vo)  -  4a(v)  (r0  -  a(v)) 


e(v)  = 


2a(v) 


When  cos  (y 0)  <  0, 

-Tocx>s(i/o)  -  \/vlcos2(yo)  -  4a(v)  (r0  -  a(v)) 


(16) 


e(v)  = 


2a  (v) 


Eccentricity  values  calculated  this  way  can  be  negative.  This  allows  us  to  plot  the  spread  of 
uncertainty  in  relation  to  the  mean  point’s  perigee  even  when  part  of  the  spread  considers  that 
point  its  apogee. 
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Finally,  we  need  an  equation  for  v(v,  tf),  which  is  the  Kepler  problem: 


[2kn  +  (E(i/) 


esinE(z/))  —  (E0  —  esinE0)] 


Where  cos  E(V) 


e  +  cos^ 

1  4-  e  cos  v 


and  E0  —  E(i/0) 


(17) 


Where  k  is  the  number  of  times  the  object  has  passed  perigee,  and  E(v)  is  the  eccentric  anomaly. 
We  write  the  angle  between  the  mean  final  position  and  the  x-axis  as  vf  which  we  write  as: 


COS  Vf 


ro 


(18) 


Keep  in  mind  that  inverse  cosine  only  outputs  values  between  0  and  n,  so  if  the  final  nominal 
position  is  below  the  x-axis,  it  needs  to  be  converted  to  equivalent  values  outside  this  range. 

A  first  order  approximation  would  be  to  assume  constant  change  of  angle  v,  which  would 
correspond  to  a  circular  orbit: 


Vi  (Y,  tf)  ~  2tt 


TP(v) 


(19) 


Where  T P  (v)  is  the  orbital  period  of  a  point  with  the  initial  position  rO  and  velocity  v.  This  leads 
us  to  write  the  following  equations: 


TP(v) 


TPq 


Q.Q 


2tt  a 
- a(v)2 

ntf  -  TPq 

2tt  4 
- an 

v7<  0 

A*r0 

VqTq  -  2 fi 


(20) 
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We  can  combine  these  with  Equation  12  to  form  one  expression  for  Equation  20: 


14 (v,  tf)  RS  27 TTltf 


/  v2r0  -  2^  \  2 

\  voro  -2 n) 


(21) 


This  approximation  will  only  hold  for  near  circular  orbits,  with  very  little  velocity  uncertainty,  and 
works  best  for  long  propagation  times,  that  is  more  than  a  third  of  an  orbit.  For  example,  for  a 
circular  mean  orbit  at  an  altitude  of  600km,  you  would  have  a  maximum  uncertainty  of  450km  at 
the  uncertainty  tails  for  a  velocity  error  of  2%  of  the  mean  velocity  after  1/3  of  an  orbit.  This  drops 
to  95  km  after  2/3  of  an  orbit.  Figure  22  shows  how  the  second  order  approximation  breaks  down 
at  the  extremes  and  the  black  line  shows  the  actual  uncertainty  cloud  after  5  orbits  for  initial  in¬ 
track  velocity  uncertainty  that  spans  from  75%  to  125%  of  the  mean  velocity  magnitude.  The  pink 
line  shows  the  equation  result  using  a  second  order  approximation  for  the  true  anomaly  in  time. 
At  the  extreme  points  of  the  uncertainty,  the  equation  no  longer  holds  true  to  the  actual  distribution. 

This  error  can  be  reduced  by  using  higher  order  approximations  such  as: 


RJ  7T  [l  —  —  e(v)^  sill  (l/L(v,  tf))  (22) 


All  together  this  will  become: 


V  (v>  tf)  ftJ  Vf  d=  [vi(v,tf)  +  l/2(v,  t/)]  (23) 

Where  positive  values  correspond  to  counterclockwise  motion  and  negative  values  correspond  to 
clockwise  motion.  For  the  same  case,  using  this  approximation  for  v(v,  tf ),  you  would  have  an 
uncertainty  of  80km  for  a  velocity  error  of  2%  of  the  mean  velocity  after  1/3  of  an  orbit  and  this 
drops  to  25km  after  2/3  of  an  orbit.  We  can  then  convert  this  from  polar  to  Cartesian  coordinates: 


For  v0-v5_  <  v  <  v0  +  vs+ 
z(v,  tf)  =  r(v,  tf)  cos  (i/(v,  tf)) 

(24) 

y(vr  tf)  =  r(v,  tf)  sin  (j/(v,  t/)) 
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Where  Equations  8,  12,  17,  23,  11,  22,  and  16  are  re-written  below  for  convenience: 


(V  tf)  -  o(v) ' 

(1  -e( 

v)2) 

1 

(  1  4-  e(v) 

*  cos (v 

CM/)) 

a(v)  = 


e(v)  = 


20-t) 


— i-qCosIVo)  ±  y/i -%cos2(i>0)  -  4a(v)  (r0  -  a(v)) 


2a(v) 


vx{\,tf)te  2 nntf 


v2r0  -  2ft 


3 

2 


Vq  r<j  - 


i/2(v,t/)  fti  I4(v)  +7T  (v/1  -  e(v)  -  l)  sill  (i/1  (v)) 


(25) 


Equation  24  is  an  approximation  of  the  in  track  uncertainty  propagation.  It  describes  the  curvature 
of  the  probability  cloud  and  also  the  in-plane  orientation  due  to  in-track  velocity  variance.  In  order 
to  use  the  equation  with  no  error,  Kepler’s  problem  would  need  to  be  solved  numerically  to  replace 
v(v,  tf). 


Uncertainty  Shape  after  5  orbits 
K  io4  75  to  1  25  velocity  magnitude  velocity  uncertainty 


Figure  22.  Equation  Extremes 
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4.0  RESULTS  AND  DISCUSSION 


Moment  estimation  is  helpful  for  Gaussian  initial  uncertainty  through  near  linear  propagations, 
but  is  a  poor  representation  of  non-linear  propagations  and  results  in  significant  information  loss, 
see  Figure  15.  Distributions  can  be  cut  for  analysis  in  an  infinite  number  of  ways,  however 
most  of  these  ways  are  not  helpful.  For  example,  to  enclose  50%  of  the  probability  of  a 
Gaussian,  the  distribution  could  be  cut  right  down  the  middle.  While  there  certainly  is  a  50-50 
chance  of  being  on  one  infinite  side,  this  is  unhelpful  in  describing  the  distribution.  We  are 
seeking  to  cut  a  figurative  banana  in  between  the  peel  and  the  good  part,  to  best  represent  the 
shape  of  the  banana  without  carrying  around  the  pesky  peel.  Another  way  to  do  this  is  with  a 
probability  threshold  slice.  This  calculates  points  that  have  a  particular  probability  in  the 
probability  distribution  function.  One  way  to  do  this  is  to  calculate  the  probability 
distribution  function  and  then  horizontally  cut  it  at  a  particular  probability  to  construct  a 
contour  of  equiprobability.  This  is  equivalent  to  finding  the  minimum  area  that  encloses  a 
particular  probability.  While  this  seems  like  a  reasonable  way  to  cut  a  distribution,  these 
properties  are  not  preserved  through  all  non-linear  propagations.  In  the  banana  picture,  this  will 
often  cut  off  the  two  ends  of  the  banana  and  leave  the  middle  part  with  the  peel.  We  propose  that  a 
standard  deviation  cut  preserves  probability  properties  through  propagation  in  full-dimensional 
states.  A  standard  deviation  cut  calculates  the  points  that  lie  a  particular  number  of  standard 
deviations  away  from  the  mean  of  the  distribution.  In  the  case  of  Gaussian  uncertainty,  this 
can  be  done  by  creating  ellipses  whose  axes  extend  to  the  corresponding  standard 
deviation.  This  type  of  cut  is  equivalent  to  a  probability  threshold  slice  for  Gaussian  probabilities, 
see  Figure  23. 

Figure  23  shows  distributions  can  be  cut  an  infinite  number  of  ways  depending  on  what  kind  of 
structure  you  are  interested  in  tracking.  The  figures  on  the  left  depict  a  Gaussian  distribution 
in  two  dimensions,  the  top  one  cut  around  one  standard  deviation,  and  the  bottom  cut  at  a 
probability  threshold  of  .08.  The  figures  on  the  right  depict  a  Gaussian  distribution  in  x  and 
a  uniform  distribution  in  y.  The  top  is  cut  around  one  standard  deviation  in  each  independent 
axis,  and  the  bottom  figure  is  cut  at  a  probability  threshold  of  .21.  The  projections  of  these  cuts  is 
shown  below  each  figure  in  red.  While  both  methods  are  equivalent  for  Gaussian  probability, 
they  are  not  for  non-Gaussian  uncertainties.  The  standard  deviation  cut  in  this  case  better 
represents  the  shape  of  the  distribution. 

We  have  developed  the  following  mathematical  criteria  for  helpful  cuts: 

1 .  The  quantitative  description  of  the  uncertainty  must  be  calculated  in  such  a  way  that  each 
uncertainty  loop  encompasses  a  particular  proportion  of  the  distribution 

2.  Any  uncertainty  calculated  in  this  way  must  be  completely  contained  by  all  uncertainty 
loops  that  encompass  a  greater  portion  of  the  distribution. 

These  can  alternatively  be  written  as: 


Approved  for  public  release;  distribution  is  unlimited. 
35 


1 .  Each  loop  as  a  function  of  probability  encloses  its  associated  probability 

2.  Each  loop  completely  encloses  all  loops  with  probabilities  less  than  itself 

We  have  found  that  with  initial  Gaussian  uncertainties,  our  method  satisfies  these  mathematical 
conditions. 


Standard  deviation  cut 


Threshold  cut 


Figure  23.  Illustration  of  Standard  Deviation  Cut  and  Threshold  Cut 
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4.1  Scenario  Testing 


Tests  have  been  conducted  on  linear,  quasilinear,  and  nonlinear  simple  scenarios  such  as  spring 
systems  and  non-linear  oscillators.  The  method  seems  to  require  no  correction  for  non-linear 
propagations  and  quasi-linear  propagations  that  only  include  error  in  a  single  dimension.  We 
verified  this  claim  by  testing  simple  linear  and  non-linear  scenarios  and  checking  it  against  a 
“truth”  Monte  Carlo  simulation.  Linear  scenarios  such  as  constant  acceleration  in  two  dimensions 
and  a  two  dimensional  spring-damper  system  were  tested.  Monte  Carlo  points  chosen  from  initial 
Gaussian  distributions  in  position  and  velocity  were  propagated  through  time.  Since  the  Gaussian 
nature  of  the  distribution  is  preserved  in  linear  systems,  we  calculated  the  mean  value  and  standard 
deviation  values  of  the  final  states  of  the  Monte  Carlo  points  in  the  x  and  y  directions  to  find  the 
“true”  standard  deviation  points.  To  test  the  hypothesis,  representative  points  were  propagated  in 
angular  increments  around  the  mean  value  for  standard  deviations  in  position  only,  velocity  only, 
and  both  position  and  velocity,  as  well  as  the  mean  value. 

There  are  an  infinite  number  of  ways  to  cut  the  distribution  such  that  a  target  cumulative 
probability  is  enclosed.  If  you  want  to  enclose  50%  of  the  probability,  you  can  do  so  by  cutting  the 
distribution  right  in  half  across  the  mean  value  and  encircling  the  remaining  distribution,  but  this 
will  not  maximize  the  information  we  gather  about  the  characteristics  of  the  entire  distribution. 
We  want  to  cut  the  distribution  as  concentric  contours  such  that  the  contour  is  preserved  through 
our  method  of  calculating  the  propagated  contour.  We  also  want  it  to  relate  to  the  probability 
distribution  function  or  cumulative  probability  function  such  that  we  can  immediately  use  a 
particular  final  point  in  space  to  calculate  collision  risk.  For  now,  we  use  a  probability  threshold 
on  the  probability  distribution  function  to  define  the  edge  contour.  Figure  24  shows  that  rather  than 
an  initial  distribution  that  follows  a  Gaussian  distribution,  the  initial  distribution  is  a  uniform 
distribution  in  two  dimensions.  The  statistical  standard  deviation  contour  converges  to  the 
calculated  standard  deviations. 


Figure  24.  Uniform  Distribution,  2-D  Mass-Spring  System 
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We  then  tested  a  quasi-linear  propagation  case,  a  propagation  that  is  linear  in  polar  coordinates 
but  non-linear  in  Cartesian  coordinates,  see  Figure  25.  The  mean-independent  equations  hold  up 
well  with  initial  tests.  Figure  25  shows  10,000  Monte  Carlo  points  were  propagated  and  their  1.5 
standard  deviation  delimitation  is  shown.  The  proposed  method  propagated  120  points,  nearly 
1/10  of  the  Monte  Carlo  Simulation,  and  give  results  that  are  smoother  and  closer  to  the  solution 
than  the  Monte  Carlo,  especially  near  the  tails,  where  the  Monte  Carlo  points  become  scarce. 
Moment  fitting  required  only  5  propagated  points  but  does  not  capture  the  distribution  well. 


First  1000  MC  points 
+  Propagated  Mean  Value 
Actual  Analytical  Delimrtaion 
Proposed  Method 
—  Monte  Carlo  Fit 
- Second  Moment  Fit 
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x  position  (units) 


Final  Position 


Figure  25.  2-D  Constant  Acceleration  in  Polar  Coordinates 

As  a  thought  experiment,  another  approach  was  experimented  with.  Each  velocity,  position,  and 
mixed  initial  points  were  propagated  for  a  range  of  alpha  values,  that  is,  for  various  standard 
deviation  distances.  These  lines  were  then  connected  and  plotted  to  form  characteristic  curves  for 
the  polar  propagation  scenario.  This  can  be  seen  in  Figure  26.  The  actual  points  for  a  particular 
direction  always  lie  on  the  mixed  lines.  Furthermore,  as  time  increased,  these  points  converged 
toV2a.  While  this  only  works  for  this  scenario  because  all  the  points  are  moving  away  from  the 
center  mean  position  and  only  contribute  to  uncertainty  in  a  particular  cardinal  direction,  the  actual 
standard  deviation  point  should  always  be  bound  by  the  velocity  only  propagation  and  the  position 
only  propagation  on  the  lower  end  of  2 a  on  the  upper  end  for  quasi-linear  and  linear 
propagations.  Since  all  the  points  are  moving  away  from  each  other,  the  position  and  velocity 
components  of  their  initial  state  become  equally  important  and  fully  mix  to  give  us  the  V2a 
value.  For  specific  applications,  this  might  be  an  easier  and  faster  way  to  delimit  the  boundaries. 
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Figure  26.  Characteristic  Curves  Polar  Projection,  2-D  Constant  Acceleration 

4.2  Probability  Risk  Assessment 


We  have  determined  that  while  this  method  is  proposed  for  better  estimation  and  assessment  of 
satellite  collision  risk,  it  can  also  be  applied  to  any  probability  risk  assessment,  especially  those 
that  require  accurate  testing  of  rare  events  and  tail  probabilities  through  non-linear  systems. 
Satellite  collision  risk  is  primarily  analyzed  in  position  and  time  space.  This  method  does  well  at 
predicting  the  rare  tail  probabilities  since  these  propagations  are  quasi-linear.  In  Section  5,  we 
demonstrate  the  mathematics  behind  our  method  is  independent  of  uncertainty  threshold  or  rarity. 
An  approximation  is  required  for  excess  folds  and  problems  that  saturate  a  space,  such  as  the  Van 
der  Pol  oscillator  shown  in  Figure  33.  For  our  analysis,  we  propagated  uncertainty  into  a  2-D 
Keplerian  propagator  and  projected  the  uncertainty  on  all  six  different  uncertainty  pairs.  The 
results  can  be  seen  in  Figure  27,  Figure  28,  and  Figure  29.  The  uncertainty  follows  the  shape  of 
the  distribution  through  folds  and  twists  that  are  not  quasi-linear. 
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In  Figure  27,  uncertainty  was  set  to  5%  of  the  magnitude  of  the  velocity  in  the  in  track  and  radial 
velocities,  and  5%  of  the  magnitude  of  the  distance  in  radial  and  in-track  position,  and  then 
propagated  via  Keplerian  two-body  motion  for  one  nominal  orbit.  10000  propagated  Monte  Carlo 
points  are  shown.  These  graphs  show  the  six  2-D  basis  loops  and  the  four  3-D  basis  loops  projected 
onto  positional  and  velocity  2-D  planes.  The  uncertainty  threshold  our  method  is  calculating  is  set 
at  .6735.  The  basis  loops  are  filling  in  the  actual  uncertainty  and  following  the  flow  of  the 
uncertainty.  These  two  projections  are  quasilinear  and  easier  to  track. 

In  Figure  28,  uncertainty  was  set  to  5%  of  the  magnitude  of  the  velocity  in  the  in  track  and  radial 
velocities,  and  5%  of  the  magnitude  of  the  distance  in  radial  and  in-track  position,  and  then 
propagated  via  Keplerian  two-body  motion  for  one  nominal  orbit.  10000  propagated  Monte  Carlo 
points  are  shown.  These  graphs  show  the  six  2-D  basis  loops  and  the  four  3-D  basis  loops  projected 
onto  the  x  and  y  phase  space  2-D  planes.  The  uncertainty  threshold  our  method  is  calculating  is 
set  at  .6735.  The  basis  loops  are  filling  in  the  actual  uncertainty  and  following  the  flow  of  the 
uncertainty.  These  two  projections  are  more  complicated  and  nodal  point  twists  can  be  seen. 

In  Figure  29,  uncertainty  was  set  to  5%  of  the  magnitude  of  the  velocity  in  the  in-track  and 
radial  velocities,  and  5%  of  the  magnitude  of  the  distance  in  radial  and  in-track  position, 
and  then  propagated  via  Keplerian  two-body  motion  for  one  nominal  orbit.  10000  propagated 
Monte  Carlo  points  are  shown.  These  graphs  show  the  six  2-D  basis  loops  and  the  four  3-D  basis 
loops  projected  onto  the  2-D  planes  labeled  above.  The  uncertainty  threshold  our  method  is 
calculating  is  set  at  .6735.  The  basis  loops  are  filling  in  the  actual  uncertainty  and  following  the 
flow  of  the  uncertainty.  These  two  projections  are  more  complicated,  as  folds  can  be  seen,  but  our 
method  is  following  the  uncertainty  carefully. 
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(a)  Position  Uncertainty  for  Keplerian  2-D  Motion 
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(b)  Velocity  Uncertainty  for  Keplrian  2-D  Motion 
Figure  27.  Keplerian  2-D  Motion 
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(a)  X  Position  Uncertainty  for  Keplerian  2-D  Motion 


Y  plot 


(b)  Y  Position  Uncertainty  for  Keplerian  2-D  Motion 
Figure  28.  Keplerian  2-D  Motion,  X  and  Y  Components 
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(a)  X  Position  vs  Y  Velocity  Uncertainty  for  Keplerian  2-D  Motion 


y  position  vs  x  velocity  plot 


(b)  Y  Position  vs  X  Velocity  Uncertainty  for  Keplerian  2-D  Motion 

Figure  29.  Keplerian  2-D  Motion,  Comparison  of  Position  and  Velocity  Uncertainties 
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4.3  Nonlinear  Dissipative  System  Analysis 


We  analyzed  simple  scenarios  including  the  Van  der  Pol  oscillator  systems  as  a  nonlinear  oscillator 
to  show  the  method’s  appropriateness  for  treating  uncertainty  in  nonlinear  systems.  The  Van  der 
Pol  oscillator  is  a  highly  non-linear  system  due  to  the  non-linear  damping  term  and  is  described 
by  the  following  differential  equation: 


9  ,  dx  k 

—  f(l  —  %2)—r  4- — x  =  0  (26) 

s  v  *  dt  m 

Where  k  is  the  spring  constant,  m  is  the  mass,  and  c  is  the  damping  constant.  When  c  is  zero,  the 
equation  reduces  to  a  simple  spring.  What  makes  this  uncertainty  distribution  so  difficult  to  track 
is  because  the  distribution  undergoes  folding  as  well  as  stretching,  much  like  the  fold-rotate-smash 
method  of  kneading  bread.  However,  when  the  distribution  loop  is  tracked  in  its  full  dimensional 
state,  these  folds  are  tracked  accurately,  albeit  with  some  resolution  error  for  longer  propagations 
due  to  the  continual  stretching.  For  these  tests,  we  used  the  Van  der  Pol  oscillator  in  two 
dimensions  with  the  following  statistics: 
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Figure  30  shows  the  full  four-dimensional  loop  projected  onto  three  dimensions  at  different 
rotation  angles.  Figure  30  graphs  show  the  three  dimensional  projection  of  a  four  dimensional 
hypersphere  propagated  through  the  Van  der  Pol  oscillator.  Views  1,  2,  3,  and  4  are  simple 
rotations  at  the  same  moment  in  time.  As  can  be  seen,  when  projected  onto  a  two  dimensional 
surface  (the  paper,  or  screen),  points  that  lie  outside  of  the  hypersphere,  can  get  projected  within 
the  two  dimensional  shape.  The  projection  becomes  the  minimum  probability  enclosed. 
Resolution  failure  can  be  seen  at  the  bottom  of  the  hypersphere  projection  in  blues  and  purples. 
In  the  axis  labels,  ‘px’  denotes  the  position  dimension  in  x,  ‘py’  denotes  position  dimension  in 
y,  and  ‘vx’  denotes  velocity  dimension  in  x. 
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We  can  also  project  these  in  all  six  combinations  of  dimension  pairs,  see  Figure  31.  Figure  31 
shows  the  uncertainty  projections  in  all  six  combinations  of  dimensional  axes  after  .5  seconds  of 
propagation,  (a)  shows  all  six  basis  loops  projected  onto  each  axes  pair,  and  (b)  shows  these 
same  basis  pairs  along  with  a  the  four  combinations  of  3-dimensional  spheres,  propagated  and 
projected  onto  each  two  dimensional  axes  pair.  The  3D  uncertainty  loops  fill  in  some  of  the  gaps 
between  the  uncertainty  basis  loops,  to  fill  out  the  edges  of  the  actual  probability  outline.  Ideally, 
propagation  would  be  limited  to  only  points  that  end  up  lying  on  this  outline,  without 
propagating  all  the  superfluous  inside  points. 


(c)  Van  der  Pol  Oscillator  View  3 


(d)  Van  der  Pol  Oscillator  View  4 


Figure  30.  Van  der  Pol  Hypersphere  Projection 
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There  is  an  interesting  scenario  when  uncertainty  in  the  phase  space  (i.e.  x  position  vs  x 
velocity),  begin  with  uncertainty  that  encloses  the  unstable  equilibrium  point  within  the  limit 
cycle  of  the  oscillator.  Mathematically,  a  point  that  begins  there  will  not  move,  and  this 
projection  method  preserves  that  possibility,  in  keeping  with  the  mathematical  criterion,  see 
Figure  32.  Figure  32  has  the  initial  uncertainty  in  the  x  position  vs.  x  velocity  projection  included 
the  unstable  equilibrium  point  (0,  0).  There  is  a  chance  that  a  point  begins  right  at  the  equilibrium 
point  and  will  stay  there  indefinitely,  and  this  method  will  keep  enclosing  this  point  within  its 
loop  indefinitely  as  well.  In  contrast,  the  y  position  vs.  y  velocity  uncertainty  did  not  include  the 
equilibrium  point. 
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x  velocity  vs  y  position 


(a)  Van  der  Pol  Basis  Loop  Projection 
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(b)  Van  der  Pol  3-D  Hypersphere  Projection 
Figure  31.  Van  der  Pol  Hypersphere  Projection  Phase  Plots 
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Figure  32.  Van  der  Pol  Enclosed  Equilibrium  Point 

Since  the  Van  der  Pol  oscillator  stretches  and  folds  the  distribution  while  keeping  it  contained 
within  a  specific  limit  cycle,  the  tails  of  the  distribution  are  also  contained  within  this  space  and 
become  folded  under  the  projection.  This  means  that  at  some  point,  the  space  will  become 
saturated  and  all  uncertainty  will  be  contained  within  any  non-zero  standard  deviation  measure, 
see  Figure  33.  Figure  33  shows  after  a  propagation  time  of  8  seconds,  the  uncertainty  is  saturated 
in  the  x  position  vs  y  velocity,  and  is  pretty  close  in  the  position  plot.  By  10  seconds,  there  is  also 
saturation  in  the  velocity  plot,  and  the  position  and  x  velocity  vs  y  position  are  very  close.  Since 
the  limit  cycle  of  a  Van  der  Pol  oscillator  contain  the  points,  the  folds  manifest  such  that  all 
uncertainty  is  contained  within  the  projection,  even  though  the  uncertainty  threshold  was  set  to 
only  1.5  standard  deviations. 
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(a)  Van  der  Pol  Projection  at  8  Seconds 
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(b)  Van  der  Pol  Projection  at  10  Seconds 
Figure  33.  Van  der  Pol  Saturated  Uncertainty 
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4.4  Consideration  of  Limit  Cycle  Solution  Flow 


Limit  cycle  analysis  allows  us  to  break  up  the  solution  flow  into  fundamental  motion  that  can  be 
used  as  non-Gaussian  basis  for  covariance  analysis.  Our  initial  analysis  consists  of  working  with 
two  non-linear  simple  scenarios  to  develop  and  test  theory  in  a  controllable  environment.  We 
worked  in  two-dimensions  tracking  x  and  y  in  position  and  velocity.  These  combinations  are:  x 
and  y  position  space,  x  and  y  velocity  space,  x  space,  y  space,  and  x  and  y  position  and  velocity 
cross  spaces  (x  position  vs.  y  velocity  and  x  velocity  vs.  y  position).  In  Figure  34  the  first 
analysis  was  a  two  dimensional  polar  to  Cartesian  coordinate  transformation.  The  Monte  Carlo 
outline  approximates  the  actual  distribution  by  enclosing  the  target  threshold  of  .67  portion  of  the 
points  in  the  least  number  of  bins  (all  other  bins,  arbitrarily  cut  to  the  same  size,  contain  fewer 
points  than  those  contained  within  the  Monte  Carlo  outline)  also  shown  in  Figure  34. 
Discrepancies  in  the  Monte  Carlo  outline  and  the  red  position  delimitation  in  the  top  left  corner 
are  due  to  only  propagating  .5  million  points.  The  red  delimitation  matches  the  actual 
mathematical  uncertainty  loop.  The  other  graphs  are  expected  to  also  have  Monte  Carlo  outline 
discrepancies  of  similar  magnitude.  Initial  conditions  were  set  as  /ur  =3.5  units,  or  =1  unit,  /j6 
=A  rad,  o9=l  rad,  /ur'  =-.5  units/sec,  of  =  1  units/sec,  ju 6 '  =.l  rad/sec,  and  06'  =.l  rad/sec. 
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Figure  34.  Polar  to  Cartesian  Coordinate  Transformation 

For  a  position  only  final  result,  the  final  distribution  does  not  depend  on  initial  velocity  spread  in 
radius  or  theta,  which  can  be  seen  as  a  collapsed  “ellipse”  in  the  position  distribution.  These 
actually  follow  the  orthogonal  basis  set  in  the  initial  polar  Gaussian  distribution  in  the  Cartesian 
coordinate  frame.  The  effects  of  these  phase  loops  in  the  other  graphs  can  be  seen  as  seemingly 
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chaotic,  but  they  follow  predictable  patterns,  which  represent  the  underlying  structure  of  the  final 
distribution. 

The  position  distribution  is  especially  easy  to  track  since  other  phase  loops  are  collapsed  in  this 
distribution.  The  red  position  only  loop  maintains  the  target  portion  of  enclosed  points  with  the 
exception  of  end  overlap:  when  an  increase  in  theta  uncertainty  is  applied,  the  ends  of  the  “banana” 
fold  over  each  other,  and  tail  probabilities  overlap  into  the  target  area  tracking  the  middle  of  the 
distribution.  This  means  that  the  red  delimitation  is  at  least  an  overestimation;  the  probability  that 
a  point  is  outside  this  area  is  equal  to  or  lower  than  the  target  proportion. 

Analyzing  Van  der  Pol  oscillation  in  the  context  of  limit  cycles  sheds  further  insight  into  the 
value  of  breaking  up  uncertainty  in  this  way,  see  Figure  35.  Figure  35  is  generated  by  setting 
spring  constants  for  x  and  y  were  set  to  1  kg/s2  each,  the  mass  was  set  to  1kg,  and  q  value  for  x 
and  y  was  also  set  to  1.  Initial  conditions  for  the  Van  der  Pol  oscillator  were  as  follows:  /ux=l 
unit,  <7X=1  unit,  qy  =4  units,  ay  =  1  unit,  and  all  velocity  initial  conditions  were  set  to  zero.  The 
same  legend  as  Figure  34  also  applies  to  these.  The  typical  Van  der  Pol  limit  cycle  can  be  clearly 
seen  in  the  center  graphs  especially  in  (c)  once  it  has  fully  formed.  This  indicates  a  loop  collapse 
that  can  be  exploited  in  this  simple  scenario.  As  time  increases,  the  distribution  becomes  more 
mixed  and  seemingly  more  chaotic,  however,  this  complexity  increases  linearly  with  time. 

We  began  by  simplifying  the  uncertainty  to  only  be  in  position.  This  allows  for  a  collapse  of  the  x 
and  y  space  into  lines,  and  to  see  the  underlying  phenomena.  While  the  previously  studied  position 
and  velocity  loops  track  stretches  in  the  distribution,  these  collapsed  mixed  space  loops  have  the 
ability  to  track  folding  of  the  distribution  over  itself.  With  these  folds  being  tracked,  the  final 
distribution  can  be  “unfolded,”  allowing  calculation  of  the  final  distribution  as  a  sum  of  its 
folded  parts.  We  do  not  expect  a  Monte  Carlo  outline  to  correctly  predict  the  structure  of  the 
distribution  in  a  way  that  can  be  tracked,  or  in  a  way  that  is  mathematically  connected  to  the 
initial  uncertainty  distribution  or  the  solution  flow  of  the  problem  for  this  oscillator,  however  a 
Monte  Carlo  outline  may  be  reconstructed  by  running  multiple  simulations  of  this  nature 
with  different  target  thresholds.  New  verification  schemes  will  be  developed  to  correctly  asses 
how  well  the  distribution  is  captured  with  these  folded  measures  in  mind.  This  will  involve 
counting  points  within  these  currently  excluded  folds  of  the  red  position  loop,  which  are 
considered  within  the  heart  of  the  distribution  by  the  orange  and  yellow  lines.  Video  analysis 
has  also  been  undertaken  to  see  the  change  in  these  folded  structures  and  to  better  understand 
how  calculations  should  properly  be  implemented. 
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(a)  Van  der  Pol  Oscillator  at  2  Seconds 


(b)  Van  der  Pol  Oscillator  at  6  Seconds 


(c)  Van  der  Pol  Oscillator  at  13  Seconds 
Figure  35.  Van  der  Pol  Oscillator  State  Space  Analysis 
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5.0  RARE  EVENT  PREDICTION  APPLIED  TO  SATELLITE 
COLLISION 


This  method  allows  for  variable  resolution  in  physical  space,  which  would  include  tail  probability 


regions.  The  mathematical  proof  shown  below  is  independent  of  uncertainty  threshold  or  rarity, 
meaning  that  tail  probabilities  are  estimated  with  the  same  accuracy  as  any  other  threshold. 
Resolution  requirements  for  tail  probabilities  may  require  additional  points,  as  these  tail  regions 
can  stretch  out  to  form  very  large  loops. 

Our  method  assumes  Gaussian  initial  uncertainties  or  any  initial  distribution  that  can  be 
transformed  to  Gaussian  through  the  use  of  a  transformation  function.  These  Gaussian 
uncertainties  are  cut  into  a  spherical  hyper  ellipse  that  are  a  particular  standard  deviation  away 
from  the  mean  point,  or  equivalently,  have  equipotential  probability  thresholds.  These  hyper 
ellipses  are  then  sent  through  the  propagation  function  and  form  post-propagation  hypersurfaces 
that  we  will  call  a  loops  or  hyper  loops.  All  hyper  loops  calculated  in  this  way  satisfy  the 
mathematical  conditions  established  above. 

We  begin  with  our  mathematical  criterion: 

1 .  Each  loop  as  a  function  of  probability  encloses  its  associated  probability. 

2.  Each  loop  completely  encloses  all  loops  with  probabilities  less  than  itself. 

Let  the  loop  function,  l(p),  be  the  enclosed  hypersurface.  Let  the  set  of  points  contained  in  lip) 
be  L{p),  i.e.  the  loop  function  is  the  boundary  of  the  points  contained  within  loop: 


lip)  =  dLip) 


(28) 


Let  the  initial  probability  density  function  be: 


(29) 


In  keeping  with  criterion  1 ,  we  define  loop  function  such  that  each  loop  as  a  function  of 
probability,  p,  encloses  p  probability: 


For  criterion  2,  we  want  a  set  of  these  loops  that  do  not  cross,  and  form  nested  sets: 


Approved  for  public  release;  distribution  is  unlimited. 
53 


(31) 


if  P1<P2,  L(pi)cL(p2) 

These  criterion  can  also  apply  for  bimodal  or  more  complicated  initial  probability  distributions, 
however  we  will  focus  on  the  case  of  Gaussian  initial  probability.  Our  method  takes  initial  nested 
hyperspheres,  ellipses  in  our  case  of  Gaussian  uncertainty,  and  propagates  these  to  full  dimensional 
hypersurfaces,  which  we  refer  to  as  loops.  Sublevel  sets  fit  criterion  1.  L(c)  is  a  sublevel  set  for / 
:  Rn  — ►  R  if: 


L(c)  =  {x  <G  R„|/(x)  <  c}  (32) 

We  would  like  to  manipulate  this  so  that  it  describes  sublevel  sets,  L(c),  in  the  final  statistics  that 
stem  from  balls,  /  (x)  c,  in  the  initial  distribution.  To  do  this  we  have  to  run  the  propagation 
backwards.  WLOG  let  the  propagation  function  be: 

f(xu  ***f*^Tt)  [2/l ■)  ***1  (33) 


If  we  index  each  final  point  to  its  corresponding  initial  point,  we  form  a  function  that  is  one  to  one: 

jf  (*^  l  j  *■*)  ;  *  *  *?  ifan]  (34) 

Since  the  function  is  one  to  one,  we  define  the  inverse  propagation  function: 

/(  -  -.!/*„)  =  [*i,-,*»]  <35) 

We  then  can  define  our  sublevel  sets  as: 


L{c)  =  {X  6  M„||  ||  /  1(yXl,...,yXn)  -  1*  ||<  c} 


(36) 


Since  our  initial  uncertainty  is  Gaussian,  c  can  be  interpreted  as  a  measure  of  standard  deviation, 
equivalent  to  probability  enclosed  for  normal  distributions: 


C  — 


(37) 
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Where  a  is  the  multidimensional  standard  deviation,  and  ap  is  a  constant  that  represents  the  number 
of  standard  deviations  required  to  enclose  p  probability: 


Xn)dX  =  p 


(38) 


This  satisfies  criterion  1.  We  expect  c  to  be  single  dimensional  inR  for  sublevel  sets;  this  can  easily 
be  done  by  scaling  each  dimension  by  its  associated  standard  deviation.  This  updated  sublevel  set 
will  be: 


L(c)  =  {X  e  M„||  ||  /  1(yXi,—,Vx  n)/a  ~  i*  II  <  Qp} 


(39) 


This  fits  the  definition  of  a  sublevel  set  for  Rn  -^R.  We  can  then  write  it  in  terms  of  propagating 
just  the  outline: 


l{p)  =  8L(p)  =  { X  e  ffin||  II  f  1(yXl t..',yXn)/o-  -  p  ||=  ap}  (40) 

Finally,  since  this  fits  the  definition  of  a  sublevel  set,  we  can  rewrite  it  forwards  as: 

l{p)  =  f(xu  ...,xn)VX||  ||  X  -  p  ||=  apa  (41) 

Which  is  the  forward  propagation  of  the  surface  points  of  a  hyper  ellipse  with  dimensions  given 
by  apa. 

There  is  an  extensive  basis  for  the  inverse  problem  in  relation  to  uncertainty  quantification 
[36],  However  its  use  in  satellite  tracking  applications  is  limited  or  non-existent.  Sensitivity 
analysis  of  the  final  distribution  to  the  mixing  of  error  in  different  dimensions  could  allow  for 
better  priority  tasking  of  ground  based  sensors.  This  combined  with  a  Bayesian  approach  to  the 
inverse  problem  could  significantly  increase  the  efficacy  of  ground  based  sensors  by  simply 
altering  tasking  algorithms  [54],  Confidence  intervals  will  quantify  the  accuracy  of 
measurements  and  calculations.  This  will  allow  mathematical  decision-making  and  help 
quantify  the  limitations  of  sensor  arrays.  Additionally,  this  will  allow  sensor  control  algorithms 
to  be  automated  at  optimal  performance. 
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6.0  CONCLUSIONS 


Understanding  how  measurements  affect  a  satellite  probability  cloud  at  a  future  time  of  interest,  and 
how  it  interacts  with  other  satellites’  probability  clouds  is  paramount  to  optimizing  sensor 
tasking  to  better  satisfy  sensor  objectives.  Accuracy  and  efficiency  could  be  increased  simply  by 
optimizing  tasking.  The  three  optimization  objectives  we  will  be  analyzing  are  localization, 
disambiguation,  and  collision  assessment.  Localization  optimization  will  require  maximizing  the 
likelihood  that  the  satellite  will  cross  the  viewing  window,  which  can  move  in  time  but  is  subject 
to  constraints.  The  problem  gets  interesting  when  multiple  satellites  are  incorporated  and  the 
objective  is  to  maximize  the  likelihood  of  seeing  all  satellites,  or  satellites  assigned  priority 
weights,  subject  to  viewing  constraints.  This  could  also  be  applied  to  multiple  viewing  windows  or 
additional  sensors.  Measurements  taken  can  also  affect  the  ability  to  localize  satellites  in  future 
viewing  windows.  Trading  off  accuracy  in  different  directions  or  dimensions  will  trade  off  the 
probability  of  localizing  some  satellites  over  others,  and  could  be  incorporated  into  observation 
algorithms  to  increase  the  satellite  siting  rate  long  term. 

Tail  probabilities  are  often  poorly  modeled,  and  as  a  result  other  methods  are  used  to  predict  rare 
events.  In  the  case  of  satellite  collision,  satellites’  mean  positions  are  analyzed  to  determine 
the  moment  of  “closest  approach”  and  a  collision  probability  calculation  is  then  made  at  that 
time  with  the  estimated  standard  deviation  of  the  positions  of  the  satellites  at  that  time.  This  is 
problematic,  because  the  “closest  approach”  time  may  not  accurately  predict  the  time  of 
maximum  probability  of  collision,  or  the  net  probability  of  collision  across  the  entire  encounter 
event.  Using  the  uncertainty  quantification  method  proposed  along  with  the  “closest  approach” 
method  of  calculating  probability  of  collision  should  by  itself  increase  the  accuracy  of  these 
predictions,  however  with  accurate  tail  probability  measures,  new  methods  of  calculating  the  time 
and  magnitude  of  the  maximum  probability  of  collision,  as  well  as  net  probability  of  collision  will 
be  possible. 

Disambiguation  optimization  would  maximize  the  difference  in  probability  of  seeing  two  or 
more  satellites.  That  is,  it  would  maximize  the  probability  of  seeing  the  target  satellite  while 
minimizing  the  probability  of  seeing  other  known  satellites  in  the  expected  viewing  space. 
When  a  satellite  is  seen,  it  often  can  be  mistaken  for  other  satellites  also  within  those  temporal  and 
localization  bounds.  We  want  to  increase  the  probability  of  correctly  identifying  this  satellite  by 
measuring  it  in  the  temporal  and  physical  spaces  that  other  satellites  have  minimal  probability 
of  being  found.  This  could  mean  adjusting  measurement  time  to  optimize  this  difference,  or 
adjust  resolution  in  position  or  velocity  along  various  directional  vectors,  subject  to 
observation  constraints.  In  addition,  maximizing  the  resolution  of  measurements  in 
particular  directions  could  minimize  the  probability  of  a  false  positive  satellite 
identification,  especially  positional  accuracy  verses  velocity  accuracy  trade-offs. 

Collision  assessment  objectives  would  require  knowledge  of  two  known  satellites  that  are  on 
potential  collision  course.  Depending  on  where  in  each  four  dimensional  (position  and 
time)  probability  cloud  each  satellites  are  expected  cross  the  other,  successive  measurements 
can  be  taken  to  prioritize  the  accuracy  of  the  probability  cloud  in  that  direction,  or  region  of 
space,  to  increase  the  accuracy  of  the  probability  calculation,  or  increase  the  chances  of 
decreasing  a  collision  probability  estimate. 
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SEZ 
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Extended  Kalman  Filter 
Probability  Distribution  Function 
South-East-Zenith 
Unscented  Kalman  Filter 
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